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I. Statement of Work 


The Air Mobility Command at Scott Air Force Base has a group of operations research 
analysts who develop and run mathematical programming models for the Command. 
Two of their most famous models are the Patient Evacuation Model and the LOGAIR 
Model. The problem generator for the Patient Evacuation Model has been placed in 
the public domain and all of the major mathematical programming software groups 
have tested their softv^'are on these models. 

Both of these Air Force models involve a large network with additional side 
constraints. For one model the side constraints describe the capacity of an aircraft 
which is shared by patients having different injury types. A burn victim may be destined 
for the burn center in San Antonio where as a soldier with a head injury may be enroute 
to the Mayo Clinic. For the other model, the side constraints enforce aircraft capacity 
for cargo sharing the same aircraft but having different origins. That is, cargo that 
originated at Tinker AFB and cargo that originated at Wright-Patterson AFB may end 
up on the same aircraft enroute to England, but this cargo must maintain its own 
identity while on this aircraft so that it will arrive at the proper destination. This is 
handled by separate networks for each origin node which are linked by mutual capacity 
constraints. In addition, the clients frequently need integer answers. 


CPLEX 3.0 has excellent capabilities, but there are many problems in this class that 
cause great difficulty for even CPLEX. The author recently ran one of these models for 
over 10 hours of cpu time on a Dec 5000/260 without obtaining a confirmed optimal 
solution. It appears, at present, that the only hope for developing robust software for 
some of these models is to exploit the underlying network structure. 


There are two manuscripts presented in this final report. The first concerns a special 
algorithm and software implementation for the constrained assignment problem. The 
second fills in a gap in the literature regarding pivot agenda algorithms when the input 
matrix is singular. Both of these papers are steps in our long term quest to solve the 
integer constrained network problem. 
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11. A Branch-and-Bound Algorithm 


The constrained assignment problem is to determine a least cost assignment of m men 
to n jobs such that an additional set of linear constraints is satisfied. This model is a 
special case of the integer network model with side constraints which in turn is a special 
case of the binary linear program. This problem is a member of the class NP- Hard and 
it is well-known that practical problems in this class are intractable. Air Force 
problems related to the assignment of pilots to schedules and the assignment of aircraft 
tail numbers to routes can be modelled as constrained assignment problems. 


The constraint matrix for the pure network problem without the side constraints has 
this wonderful property of being totally unimodular. Hence, every basis has 
determinate equal +1 or —1, every basis is triangular, and every extreme point has 
integer components. This problem is a member of the class P and is relatively easy to 
solve. By appending only a single side constraint, the unimodularity property is lost, 
bases are not triangular, extreme points may not be integer, and the problem is a 
member of the class NP-Hard. Unfortunately, almost all real-world problems have 
one or more side constraints. For Air Force models, it is common to have some type of 
aircraft capacity or budget constraint. 


The objective of this study was to develop and empirically evaluate a new algorithm for 
this model. The algorithm relies on the branch-and-bound strategy and is designed 
for problems having a large network and a limited number of side—constraints. It 
exploits an algorithm that we developed earlier for the assignment problem having a 
single side constraint. This work uses an excellent assignment code that our research 
team developed and handles the side constraint via the use of Lagrangean relaxation. 


The paper resulting from this study appears in Appendix A of this report. It has been 
submitted for publication and is currently under review. 



III. Recovery from Numerical Instability 


In my opinion, the best technique available for solving constrained networks is to use a 
simplex based algorithm in which the basis is partitioned into two parts, one part 
associated with the network constraints, and one part associated with the side 
constraints. The component associated with the network part can be maintained as a 
rooted spanning tree and all operations involving the inverse of this component can be 
executed using specialized labeling algorithms. Another component, corresponding to 
the side constraints is called the working basis. It is the inverse of this working matrix 
which is needed for the operations required by the simplex method. 


In our system, the inverse of this working basis is maintained in factored form and every 
pivot involves the addition of either one or two new factors to the eta file. Periodically, 
say every 50 to 100 pivots, the working basis is reinverted and a new eta file is developed. 
The new eta file is smaller than the old one and much of the round-off error which has 
been introduced during the previous pivots will have been eliminated. 

The first step in the procedure to obtain a new factorization is to determine a 
permutation of the rows and columns so that the sparsity property of this matrix will be 
maintained in the factorization of its inverse. In the literature, the permutation of the 
basis is known as the pivot agenda and there are several algorithms for obtaining a good 
pivot agenda. All pivot agenda algorithms assume that the input matrix is nonsingular. 

In our work with specialized partitioning methods for networks with side constraints, 
we discovered that due to the nature of the side constraints a singular input matrbt 
would eventually be presented to the pivot agenda algorithm. When this occurs, all of 
the pivot agenda algorithms, of which we are aware, fail. The objective of our 
investigation was to present recovery procedures, using a variation of the 
Hellerman - Rarick P3 algorithm, for the case in which the input matrbc is singular. The 
results of our investigation are presented in Appendix B of this document which has 
been submitted for publication and is currently under review. 
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ABSTRACT 


This manuscript presents a branch-and-bound algorithm to obtain a near 
optimal solution for the constrained assignment problem in which there are only a 
few side constraints. At each node of the branch-and-bound tree a lower bound is 
obtained by solving a singly constrained assignment problem. If needed, Lagran- 
gean relaxation theory is applied in an attempt to improve this lower bound. A 
specialized branching rule is developed which exploits the requirement that every 
man be assigned to some job. A software implementation of the algorithm has 
been tested on problems with five side constraints and up to 75,000 binary vari¬ 
ables. Solutions guaranteed to be within 10% of an optimum were obtained for 
these 75,000 variable problems in from two to twenty minutes of CPU time on a 
Dec Alpha workstation. We believe that this is the current best algorithm and 
software implementation for the constrained assignment problem having a limited 
number of side constraints. The behavior of the branch-and-bound algorithm for 
various problem characteristics was also studied. This included the tightness of the 
side constraints, the stopping criteria, and the effect when the problems are unbal¬ 
anced having more Jobs than men. 
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I. INTRODUCTION 

The constrained assignment problem is to determine a least cost assignment 
of m men to n jobs such that an additional set of constraints is satisfied. This 
model is a binary linear program and may be stated mathematically as follows: 




i 


» 


» 




minimize Y CyXy 


(1) 

(i.j) 6 A 

subject to ^ Xjj 

= 1 , i = 1, ...,m 

(2) 

j; (i,j) e A 

Z Xij 

< 1 , j = 1, ...,n 

(3) 

i : (i. j) S A 

Xij 

€ {0.1}, all (i,j) € A 

(4) 


< r*', k=l,...,s 

(5) 


(i.j) e A 


where xy = 1 implies that man i is assigned to job j at cost of Cy, d-j denotes the 

coefficient of Xy in the kth side constraint, r’^ denotes the right-hand-side for the 
kth side constraint, and A is the set corresponding to the feasible assignments. 
Note that the problem allows for more jobs than men. Many practical problems 
have this feature. It also allows for the case in which the number of men exceeds 
the number of jobs. For this case, one simply reverses the definition of men and 
jobs. 
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Since (l)-(5) is a binary linear program, all the literature on integer pro¬ 
gramming applies (see Geoffrion and Marsten [5], Salkin [12], Parker and Rardin 
[11], Nemhauser and Wolsey [10]). In practice most integer programming models 
are either solved as a linear program and the solutions are rounded using some 
heuristic or branch-and-bound is used in an attempt to obtain a solution within a 
prespecified tolerance. 

A special case in which the side constraints have the generalized upper bound 
(GUB) structure has been studied by Ali, Kennington, and Liang [2]. A relaxation/ 
decomposition procedure that involves solving a series of pure assignment prob¬ 
lems is used successfully. Ball, Derigs, Hilbrand, and Metz [3] also present an 
algorithm for the matching problem with generalized upper bound side constraints. 

Another special case for s=l and m=n has been studied by Gupta and Sharma 
[6], Aggarwal [1], Mazzola and Neebe [9], Bryson [4], Kennington and Moham- 
madi [7,8]. The only specialized algorithm for (l)-(5) is the two phase procedure 
of Mazzola and Neebe [9]. The first phase uses subgradient optimization to obtain 
an advanced start for the branch-and-bound method used in the second phase. 

The objective of this study was to develop a new branch-and-bound algo¬ 
rithm to solve the constrained assignment problem and to provide an empirical 
analysis of this algorithm on a variety of assignment problems having only a few 
side constraints. All empirical analysis was performed on problems having five 
side constraints. 
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n. i HE BRANCH-AND-BOUND ALGORITHM 
The branch-and-bound method can be viewed as a divide and conquer strat¬ 
egy. If a problem cannot be solved, then it is partitioned into several smaller prob¬ 
lems. The best of the solutions to the smaller problems will be the solution to the 
original problem. 

Consider the problem P(S) s min{ cx: x g S}. Using the terminology of Geof- 
frion and Marsten [4] let v[P(S)] denote the optimal objective function value for 
the problem P(S). Let x* denote an incumbent for P(S) with objective value of v'. 
Let F(S) denote a relaxation of P(S) and let CL denote the candidate list. The 
generic branch-and-bound algorithm may be stated as follows: 

Input: 

1. The problem, P(S). 

Output: 

1. The solution vector, x'. 

2. The objective value corresponding to x*, v’.(v* = oo implies that S = O.) 

Procedure BAB; 

Begin 

initialize: 

CL:= {P(S)}, V:= «; 
while CL 4> do 

comment: select a candidate problem for analysis, 
select P(U) € CL, CL:= CL \ {P(U)}; 
if F(U) has a feasible solution, then 
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if v(P(U)] < v‘, then 
let X be an optimum for P(U); 
if X G S, then 
x‘:= X , v';= cx; 
else 

apply a heuristic to x in an attempt to create x such that x g S and 

cx < v‘; 

if successful, then x;= x, v':= cx; 
comment: branching 

create U,, Uj, Up such that UiuUau ...uUp = U and U,nU, = ‘J’ 

for all 15“] G {1. 2.p}; 

CL;= CLu{P(U,). P(U 2 ). .... P(Up)}; 

end if 
end if 
end if 
end while 
end. 
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ffl. AN EXACT ALGORITHM FOR 'ITiE CONSTRAINED 
ASSIGNMENT PROBLEM 

In this section we present a specialized branch-and-bound algorithm for the 
constrained assignment problem. The constrained assignment problem is a binary 
problem there.^ore at each node of the branch-and-bound tree either an assign¬ 
ment is prohibited by fixing the corresponding variable to zero or a man is perma¬ 
nently assigned to a job by fixing the corresponding variable to one. 

3.1 The Relaxation 

Let D denote the matrix corresponding to the coefficients in (5), x denote the 
vector corresponding to the binary decision variables, c denote the vector of costs, 
and r denote the vector of right-hand-side values for the side constraints. Then the 
constrained assignment problem may be denoted as P(S) where S = {x: (2), (3), 
(4), (5)}. Let 1 denote a vector of I’s and S = {x: (2), (3), (4), iDx < Ir}. Then 
P(5) is a valid relaxation for P(S). Application of the algorithm in [7] to solve the 
singly constrained assignment problem will yeild a lower bound for P(S), the opti¬ 
mal Lagrangean multiplier corresponding to the constraint IDx < ir, and x g S. If 
JT G S, then cx is an upper bound for P(S). 

Let S = {x; (2), (3), (4)}. Recall that a Lagrangean dual for P(S) is the 
problem max {L(a) :a > 0} where L(a) = min {cx+a (Dx-r): x g S }. We use the 
optimal Lagrangean multiplier and x from the singly constrained algorithm to 
form an advanced starting value for a. Let y denote the solution for L(a). Then 
2 = Dy-r is used to modify a for successive steps. A limited number of these steps 
will be performed in this algorithm. 
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3.2 The Branching Rule 
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Consider any node in the branch-and-bound tree. If the relaxation, P(S) has 
no feasible solution, then this node may be fathomed. Otherwise, an assignment 
will be used to create the branches as illustrated in Figure 1 . 

Figure 1 here 

For a given node in the branch-and-bound tree let F‘= {(i,j): x^ = 1 } and F°= 
{(i,j): Xjj = 0}. Let U= {x e S: xy = 1 for all (i,j) e F‘ and Xy = 0 for all (i,j) € P}. 
The relaxation solved at each node in the branch-and-bound tree is P(U) which is 
a singly constrained assignment problem with some assignments fixed. Let x g U , 
T = {(i,j) : Xij = 1 , (i,j) g F‘}, and t = iTj. Consider the t +1 subsets of U (U,, U 2 , 
...,i7,.i) created in the following manner: 

Ui = { X G tJ: Xijjj = 0, (i„j,) G T}, T, = T\{(ii,j,)}; 

U 2 = { X G tJ: Xi,jj = 1, Xiyj = 0, (i 2 ,j 2 ) G T,}, T 2 = T,\{(i 2 ,j 2 )}: 

tJa = { X G U; Xi,j, = 1, Xyj = 1, Xi,j3= 0, (ia.ja) G T 2 }, T 3 = T2\{(i3,j3)}; 

U, = { X G U: Xi,j, = 1 , ..., Xi,.,j^_j =1, Xij, = 0, (i,,j,) G T,.,}; 

U,*, = { X G U; Xjj = 1 for all (i,j) g T} . Note that, U = U,uU 2 U ... uU,*, and 
U’intJj= <P for all i j. 

Consider a node in the branch-and-bound tree having f, edges fixed at 1. 
Then branching from this node will produce t+1 = n-L+l new candidate problems. 
From Figure 1 it can be seen that the last node (i.e., U,»,) need not be created 
since it was examined at the parent node. Therefore, each branching produces t 
candidate problems. 


I 
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3.3 The Candidate List 

For our implementation of the branch-and-bound algorithm, ASSIGN+1 [7] 
will be applied to P(tJi) and the results placed in the candidate list (i.e., a problem 
is solved before it is placed in the candidate list). The motivation for placing solved 
problem in the candidate list is that the solution for P(i7i) can be easily modified to 
obtain an advanced starting solution for P(Uj*i). Therefore solving the sequence of 
problems P(Ui), P(U 2 ), ..., P(Ut) should require only a moderate amount of com¬ 
putational effort. Hence, each entry in the CL consists of the five tuple (P, FJ x , 
P, u) where x e Ui, ^ < v[P(Uj)], and u the optimal Lagrangean dual for the 
singly constrained problem. 

3.4 The Fathoming Rules 

At any node p of the branch-and-bound tree let U, F', and F“ be as defined 
in Section 3.2. Let M = {1, 2, 3,..., m}\ {i: (i,j) e F^}, then the following rules may 
be used to fathom a node. If 

Z Z r”in(dl) : (i, j) e A \ F°) > r'^ 

(i, j) e F’ i € M 

for any k, then node p can be fathomed. That is, no selection of the free variables 
will satisfy the kth side constraint. If 

2] Cij+ Z rn>n(Cij : (i.j) G A \ P) > v’ 

(i,j) 6 F* i e M 

then node p can be fathomed. That is, no selection of the free variables will result 
in a solution superior to the incumbent. 

If min {IDx: x e U } > Ir, then node p can be fathomed. That is, no selection 
of the free variables will satisfy all side constraints simultaneously. Let ^ be the 


A-9 






best lower bound obtained for node p and e be the termination tolerance. We will 
fathom node p if v‘ -/S < c/S. Using this rule with € = 0.1 results in a solution from 
the branch-and-bound procedure guaranteed to be within 10% of the optimum and 
€ = 0.01 produces a solution within 1% of an optimum. 

3.5 The Algorithm 

In this section we combine the information presented in Sections 3.1-3.4 to 
construct the ASSIGN+s algorithm. 

Input: 

1. The cost vector, c. 

2 . The feasible region S. 

3. The set of (man, job) pairs corresponding to eligible assignments, A. 

4. Termination tolerance, e. 

5. The maximum execution time, tmax. 

6 . The maximum number of Lagrangean relaxations to be solved at each node, 
limit. 

Output: 

1 . The solution vector, x'. 

2. The objective value corresponding to x’, v‘. (v' = oo implies that the problem is 
infeasible.) 

Procedure ASSIGN+s; 

Begin 

initialize: 

comment: Node 1 in the Branch-and-Bound tree. 
v-:= oo, P:= <P, F‘:= <t; 
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ASSIGNP1(P(S), X, /S, u); 

if /? = -oo then terminate; 

if X s S then x':= x , v’:= cx; 

else LAGRANGE(5c , u); 

comment: Tolerance test for fathoming. 

if v‘ - ^ < €/S then terminate; 

CL:= {(P,F, X. u)}; 
while CL 5 ^ do 

SELECT A PROBLEM(F°.P. x, /S. u); 

BRANCH(CL, F“,F‘, x, u); 

end while 
end. 

procedure ASSIGNP1(P(D), x. /?, u); 

Begin 

initialize: 

-00; 

apply the ASSIGN+1 algorithm (Kennington, Mohammadi [1991]) to P(U); 
if P(t7) has a feasible solution then 
let X be the best feasible solution found for P(U); 
let P be the best lower bound foimd for P(U); 

let u be the optimal Lagrangean dual for the singly constrained problem; 
end if 
end. 
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procedure SELECT A PROBLEM (P,F^ x, jS, u); 

Begin 

select (P.F, x, yS, u) e CL, CL:= CL \ (P,F\ x, u); 

end. 

procedure BRANCH(CL, P.F^, 5c, yS, u); 

Begin 

initialize: 

t;= 0, G;= = 1, (i,j) g F>}; 

M:= {1, 2. 3, ..., m}\ {i; (i,j) e F}; 
while G ’“O do 

comment: Fix a variable at zero, 
let (i„j,) e G, G:= G \ P:= P U 

ASSIGNPl(P(iJ), 5c, /3, u); 
if yS * ~<x) then 

if X e S and cx < v* then x';= x , v':= c5c; 

else LAGRANGE(5c , yS, u); 

comment: Tolerance test for fathoming. 

if V - y3 > eyS then CL;= CL u {(F°,F\ 5c, yS, u)}; 

comment: Permanently assign man i, to job j,. 

M:= M N {i,}, F‘:= F‘ u {(ii.j,)}; 

P:= P u : (ii,j) € A for all j} u {(i,ji) : (i,ji) € A for all i }\ {(ii,j,)}; 
comment: Assignment polytope feasibility tests. 

^ X X : (ij) € A) > v’ then return; 

Ojyepi iEM 


A-12 





i 


m 


for k=l,...,s 

if ^ X • 0> j) ^ return; 

(i,j)eF‘ iGM 

end for 

end if 

end while 

end. 

Procedure LAGRANGE(x , y3, u); 

Begin 

initialize: 

a:= ul, y;= x , t:= 1; 
while (v*-yS > and t < limit) do 
z:= Dy-r; 
for i=l,,,.,s 

if Zk > 0, then ak.= 1.25a*; 
if Z/c < 0, then a*;= 0.75a*; 

end for 

let y be a solution for L(a) and v[L(a)]}; 

if y e S and cy< v' then x':= y, v*:= cy; 
t;= t+1; 
end while 

end. 
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This branch-and-bound algorithm exploits the structure of the model (l)-(5) 
in several ways. Permanent assignment of a man to a job implies that all other 
variables involving this man and job may be fixed at zero. The relaxation was a 
singly constrained assignment problem for which near optimal integer solutions 
can be obtained using the results in Kennington and Mohammadi [8]. Special fath¬ 
oming rules were developed which were based upon the assignment polytope. 



IV. EMPIRICAL ANALYSIS 


► 
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The algorithm ASSIGN+s has been implemented in software and empirically 
analyzed on an Alpha workstation by Digital Equipment Corporation. The code is 
written in Fortran and uses ASSIGN+1 (see Kennington and Mohammadi [7]) to 
solve the singly constrained assignment problems. ASSIGN+1 is an implementation 
of the Lagrangean relaxation algorithm for sparse singly constrained assignment 
problems. 

We developed a test problem generator with the following inputs: (i) the 

number of men, (ii) number of jobs for each man, (iii) the maximum cost, c, (iv) 
the number of side constraints, s, and (v) the side constraint multiplier, k. Both the 
costs and the side constraint coefficients are uniformly distributed over the range 

[0, c). We randomly generate a feasible assignment, K, and set the right-hand- 

side of the side constraints, r, to kDlc. Obviously, as k becomes smaller, the feasi¬ 
ble region becomes smaller and for sufficiently small k {x: (2), (3), (4), and (5)} is 
usually empty. 

The generator was used to generate two sets of 400x400 problems described 
in Table 1. As Table 1 indicates, the problems generally become more difficult as 
k becomes smaller and for k small enough the feasible region is empty. For all 
runs, the stopping criteria used is €=10% and the % deviation reported in column 8 
gives a guarantee on the deviation from optimality. All times are the CPU time and 
exclude the time for both input and ouq)ut. The run with problem number 2 having 
k=0.4 was terminated after the candidate list grew to 25,000 entries. As we ex¬ 
pected, there exists problems which cannot be solved in a reasonable amount of 
time and storage using this approach. Tightly constrciined problems having 48,000 
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binary variables definitely stretches the capability of this software implementation 
of our branch-and-bound algorithm. 

Table 1 here 

Tables 2 and 3 give our empirical results with 30 randomly generated assign¬ 
ment problems with various sizes all having five side constraints. For all of the 
problems tested we were able to find a solution guaranteed to be within 10% of an 
optimal solution. The six smallest problems have 3,000 binary variables. Five of 
these were solved in less than two minutes each and one required about six and 
one-half minutes. The six largest problems had 75,000 binary variables and were 
all solved in less than twenty one minutes each. The most difficult problems 
(300x300) have 27,000 binary variables. Two of these six problems required eighty 
minutes to solve. These two difficult problems also had very tight side constraints. 

Tables 2 and 3 here 

This work was motivated by models for assigning sailors to ships and for this 
application the number of jobs always exceeds the number of sailors available. 
Frequently the job list covers a longer period than the list of available sailors which 
produces a large imbalance in n and m. Tables 4 and 5 present our empirical 
results from solving 18 unbalanced assignment problems with five side constraints. 
For the 300x600 problem with k=0.6 presented in Table 5, the run was terminated 
due to candidate size limit. For all other test problems we were able to obtain a 
solution within 10% of an optimum. 

Tables 4 and 5 here 

For all test problems we search for a solution within 10% of an optima. To 
study the effect of the tolerance value on the performance of the algorithm we 
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solved two 200x200 and two 200x400 problems with different tolerance values. 
Figure 2 indicates that, as expected, a decrease in the tolerance value leads to an 
increase in the execution time. For all four problems, a point was reached in which 
a slight decrease in the tolerance resulted in a large increase in the solution time. 

Figure 2 here 
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V. SUMMARY AND CONCLUSIONS 

We have presented a branch-and-bound algorithm for the constrained as¬ 
signment problem. The algorithm is applicable for both balanced and unbalanced 
assignment problems having inequality side constraints. The algorithm uses a spe¬ 
cialized branching rule that exploits the underlying structure of the problem. 
Bounds are obtained by solving a singly constrained assignment problem followed 
by a few iterations with a Lagrangean relaxation. 

We presented empirical results for both balanced and unbalanced problems 
having five side constraints. For problems having 75,000 binary variables, solu¬ 
tions guaranteed to be within 10% of an optima were obtained in less than twenty 
one minutes on a Dec Alpha workstation. Our analysis indicated that as the side 
constraints become tighter the execution time and number of branch-and-bound 
nodes increases. For one of the 300x300 problems having 27,000 arcs, the execu¬ 
tion time increased from about one minute to about eighty minutes as a result of 
side constraint tightening. Our analysis also indicates that the performance of the 
algorithm on unbalanced problems is generally better than its performance for the 
balanced problems with the same number of binary variables. The Navy personnel 
assignment problems which motivated this study are all unbalanced models. 

For problems of this type, having only a few side constraints, we believe that 
this is the current best algorithm and software implementation available. Solutions 
guaranteed to be within 10% of an optimum should be obtained for most problems 
having fewer than five side constraints and fewer than 20,000 arcs. However, as 
with any other branch-and-bound based procedure, we found difficult problems 
which required an extraordinary amount of computer time. 




VI. REFERENCES 


1. V. Aggarwal, “A Lagrangean-Relaxation Method for the Constrained 
Assignment Problem,” Computers and Operations Research vol. 12 pp. 97-106, 
1985. 

2. I. Ali, J. Kennington, and T. Liang, “Assignment with En Route Training of 
Navy Personnel," Naval Research Logistics Quarterly vol. 40 pp. 581-592, 

1993. 

3. M. Ball, U. Derigs, C. Hilbrand, and A. Metz, “Matching Problems with 
Generalized Upper Bound Side Constraints,” Networks vol. 20 pp. 703-721, 
1990. 

4. N. Bryson, “Parametric Programming and Lagrangian Relaxation: The Case 
of the Network Problem with a Single Side-Constraint,” Computers and Opera¬ 
tions Research vol. 18 pp. 129-140, 1991. 

5. A. Geoffrion and R. Marsten, “Integer Programming Algorithms: A 
Framework and State-Of-The-Art Survey,” Management Science vol. 18 pp. 
465-491, 1972. 

6 . A. Gupta and J. Sharma, “Tree Search Method for Optimal Core 

Management of Pressurised Water Reactors,” Computers and Operations Re¬ 
search vol. 8 pp. 263-269, 1981. 

7. J. Kennington and F. Mohammadi, “The Singly Constrained Assignment 

Problem: An AP Basis Approach,” Technical Report 93-CSE-25, Depart¬ 
ment of Computer Science and Engineering, Southern Methodist University, 
Dallas, TX 75275, 1993. 





8. J. Kennington and F. Mohammadi, “The Singly Constrained Assignment 

Problem: A Lagrangean Relaxation Approach,” to appear in Computational 
Optimization and Applications. 

9. J. Mazzola and A. Neebe, “Resource Constrained Assignment Scheduling,” 

Operations Research vol. 34 pp. 560-572, 1986. 

10. G. Nemhauser and L. Wolsey, Integer and Combinatorial Optimization, John 
Wiley and Sons: New York, NY, 1988. 

11. R. Parker and R. Rardin, Discrete Optimization, Academic Press Incorporated: 
New York, NY, 1988. 

12. H. Salkin, Integer Programming, Addison-wesley Publishing Company: 


Reading, Massachusetts, 1974. 










A-21 








lime 



a. Problem #1 (200x200) 
with 1200 arcs 


Time (min.) 


b. Problem #2 (200x200) 
with 1200 arcs 


Time (min.) 




c. Problem #3 (200x400) 
with 1200 arcs 


d. Problem H (200x400) 
with 1200 arcs 


Figure 2. Plots of time versus the stopping tolerance for four problems. 
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Table 1. Empirical results from branch-and-bound algorithm for 400x400 assignment 
problems with five side constraints (48,000 columns and 805 rows). 


Prob. 

B 

# BAB 
Nodes 

# AP’s 
Solved 

Time 

(min.) 

LB 

UB 

% 

Deviation 

Node # of 
Incumbent 


1.0 

189 

216 

0.71 

5,270 

5,309 

0.74 

13 


0.9 

244 

1,958 

4.74 

5,505 

5,768 

4.78 

22 


0.8 

274 

2,119 

3.96 

6,999 

7,271 

3.89 

78 

1 

0.7 

952 

8,288 

12.66 

10,991 

11,725 

6.68 

859 


0.6 

14,142 

124,222 

166.00 

20,820 

22,713 

9.09 

13,943 


0.5 

4,214 

35,362 

45.39 

47,446 

50,343 

6.11 

4.157 


0.4 

1 

2 

0.00 

problem has no feasible solution 


1.0 

224 

503 

1.66 

5,370 

5,559 

3.52 

2 


0.9 

253 

1,801 

4.57 

5,830 

6,132 

5.23 

34 


0.8 

700 

5,434 

9.26 

7,815 

8,373 

7.13 

492 


0.7 

3,114 

25,786 

42.93 

12,510 

13,649 

9.10 

2,938 

2 

0.6 

3,094 

28,694 

40.92 

23,606 

25,188 

6.70 

2,871 


0.5 

3,505 

29,642 

31.45 

50,939 

54,690 

7.36 

3,340 


0.4 

25,000^ 

79,160 

53.37 

148,825 

no feasible solution obtained 


0.3 

1 

2 

0.0 

problem has no feasible solution 
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Table 3. Empirical results with problem set 2 using the branch-and-bound algorithm. 

(All problems have five side constraints.) 



A-25 































































































































































































































































T('{ liiii( al Rcixtrt 94-C'SE-32 


Recovery from Numerical Instability 
During Basis Reinversion 


Jeffery L. Kennington 
and 

Riad A. K. Mohamed 


DoparTinoiit of Coiiipiucr Science and Engineering^ 
School of Engineering and Apjjlied Science 
Sonthern Methodist Elli^■ersitv. 

Dallas. Texas 7 j275 


August 1994. 


B-l 











Abstract 

All ()! 1 iic pi'caxiuiK'il |)i\(.; in-ciicla aluutit Inns t ha I extend i 1 k' llellcrmai/- 
Hai ich I’d aluoi it Imi assuna' i liai i lie iiipiii inairix is iionsingular. Due to ini- 
niei ii ai ill'! allilii i hi' a"iini|ii ion may lx- x'iolated and i lie'i' aliiui ii lims iail. 
We jireseni a moditn atidii ut the I'd aluorilhm which inchides a jncxeduic 
t<< recocei IVfmi ihi- type of niiiwriiii) insial>i)il,\', ] lie )-ec(j\'ei'\' piuci'dnre i' 
inteiirateil into Pd in 'uch a wav tha: all imndoiis work can he maintained 
an<l it reduce- the likeliheod that addiiioiial recovery will he retjuired. 
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1. Introduction 


III liiK'rti [)ixj,£>irtiiiniiiii; "X'tnii'. iiii imj»urt^iii rominMH’iii i' l lie hlnoin hm tm iil> 
t ainiiiii « iicw liii turi/,ii ioi, toi t In- itivor-c »>!' i he J lie iii>i 'i r]) in i li<' aliiui il inn 

i> til iiciniui at iiiii oi 1 tic luws and <()luinn> ul the lia-'i'-. 'u that i lie '•jia; 

^i1y )jioi)criy oj the hasi^ will l>c inaiiitaincd in tin lai'torizai mn ul it > in\ci>c. In 
the lit nature, the iienmit at ion of the ha^ii^ is known as a pivot lup iidti. and there are 
several algorithms for ohtaining a good pivot agenda. 

In 1III i. I lelleinian and Hariek '> int rod need t he preas>ii;iied pi \ ui pioreiliiri ■ ' 1 ’i-)' 
ior olitainiiiii '-inh a iiermutatiun. In lUT'd. th<‘\‘ added an initial •'urt that jieiinuie' 
ill' mat ri:-. I" iuwei hiuel-; t riaii'inla: lurm. a in! then applies a 'implihi'd \ ei sii ,n ,i n, 
Pd al'.;>iril hill lu each liioel; 1 he\ ea'ie'l thi' ahjurillini the jiart it iuiiei] pii'assi!rneil 
|)i\oi procedure iPli. see (I. Dull .2. presented a siiipile altiorithiii lor periiiuiini; 
the matrix so that its diagonal ha' a minimum niimlx'r of zeros. I'.risnian it <il. d . 
disnis'<-d the po"ihilit v of a't ruet ii rally zero pivot arising in Pi. ainl sug'^e^i ed a vari¬ 
ant tiiat avoid' thi' prohlem. called the piecautionary |jartiticuied juias'i'jiied pivoi 
procedure ip.'i. .Xrioli (t al. 1 . re|)orted that the P." al'iorilliin doc' not a'hlre" 
the luohleiii' a'sociati'il with 'iiiaii pivots, ainl thai Pi jrerlormed heller than Pd 
when takiii'i intoaciount numerical j^ivotinti. Halleisley and .Macklev 1 . desnihed 
a transposed version of Pd which permutes the matrix into an up];ei triangular form 
with row spikes extendiiie below the diaeonal. and sugeested a technique whicli alter¬ 
nates between the classic and transposed versions of P3 for reducing the build-up of 
nonzeros during factorization. Sankaran [7j. presented some new results on optimal 
spike configuration, using an ordering heuristic for doubleton columns. 

.\I1 of these methods assume that the input matrix has at least one nonzero entry 
in every row. In our work with specialized ))artitioning methods for networks with sid<“ 
constraints, one maintains the inverse oi a working Ijasis. Q. corresjionding to the side 
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cuiii't limits. 1 )u<' t () t lie iiai UK' ui r»'rtaui i vi«'^ i il side const laint ^ <mi 1 liiuti • | ui i isn>11 
a CJ can iic (ic\cliijicd iiavinc one or tnoic rows, i-acli entrv of wliii h k sniallt i liiai' 
our zero toltuancc. Ht'iicc. tlic usual pivot agenda algoritlmis fail. For prolilenis in 
this class, we discovered that nuinerical instabilities frequently occur. Scaling can 
alleviate some of these problems, but can not guaratitee that this difiiculty will not 
arise. \\ hen a row of all zeros (or near zeros) is discovererl. one generall.v K'places some 
column of Q with an artificial column and the simplex method is continued. This 
idea is referred to as recovery (i.e. the algorithm recovers frotii numerical instability). 

The objecti\'e of this investieation is to nresent reccjveiy procedures, usinc a \ari- 
atit of the Hellerman-Karick !'.{ algorithm, for tin* case in whi<'h the input matrux 
has vine or more rows ol all zr'tos. Our variation allows decomjrosing the bum]> into 
smaller blocks, and permits interclianses among the una.ssigned rows to avoid what 
Erisman ft al. refer to as a .-^tnicturalh) z>r<> pix'ot. Recovery involves rejrlacing some 
column of the input matrix with an aj>))ro))riate column of the identit.v matrix. Tlu' 
difiii'iilt part is tr* determitie the column to be reirlaced so that the current work on 
the |)ivo1 agenda can be retained. .An arbitrar.v selection of the column to be- reiilaced 
can lead to failure of the f‘d algorithms. This failure occurs when it is discovered that 
there exists a row having all zeros in unassigned columns, or a structurally zero juvot 
may arise. Both of these conditions will l>e demonstrated in the examjile to follow. 

Consirler the matrix illustrated in Figure la. which has all zeros in row 10. Re¬ 
placing column 3 by elO ( a column with a nonzero in row 10 and zeros elsewhere) 
and applying the P3 algorithm results in all zeros for unassigned columns in row 3. 
as illustrated in Figure lb. This gave us the intuition to confine the column replace¬ 
ment to the lower triangular blocks. However, this does not completely solve the 
problem. Consider the matrix illustrated in Figure Ic. which also has all zeros in 
row 10. Replacing column o by elO. then applying the P3 algorithm, will introduce 
a structurally zero ])ivot at location (T.S). as illustrated in Figure Id. 









n 


[ 



! 


X 

X X 
XXX 
X XX 
X XX 

X XX X 

X XXX 

X X X X X 

xxxxxxxxxx 



A new zero 

row 


a aatrlx as parautsd by P3. 
column 10 is auggsstad for rsplacmanc. 



(b) whan raplacing column 3 by alC. 



(c) a block as permutad by P3, (d) whan raplacing column 5 by alC. 

column 10 Is suggastad for raplacmant. 


Figure 1: I-xaiujile^ of failui'' (iuo tu iucunfci coluiiiii !<’i)la( (Miii')ii. 

2. A Pivot Agenda 

L<’l B Ix' a l^ooleaii maliix oi ludci n. aixl lei h, deiioU' llic ij" ('IciikmiI id B. 
A ]>i\oi aseixla for B i*- a iieiniulatiuii of ilic rows and rolnniii'- ol B whii li vii'hK a 
inatrix lla^■in'l soriH’ dcsiralrlc j^rojxTty. For tiiis work the desiiahle jjrojjeriy i> lowei 
triangular form. If this is impossible, then we seek a near lower triangular form. Let 
C and 7? denote the sequence of v columns and rows in a pivot agenda. That is. the 
column (row) in the permuted matrix is CaITv*-). 

.A pivot agenda can be developed lyv using three distinct procedures. One proci'- 
dure searchs for column singletons and places the.se columns at tlie end of the pivot 
agenda. .Anothei procedure search^ for row singletons and place> these row- at the b< - 
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uiiiiiiiit; (<1 1 1 / 1 ' 1 li'' t liii 'l |»r>M cihiic ilctcriiiiiic- i Im' |)|'i iniii .ii in!; In: i i;' 



rciuaiiiiiiii inv.' aiii! < niuiiiu^ to < on-^t I'lict one <ir nmn' lowci liiaiititiiai l/lix I;-. I ti<-, 
j)rof<'(lui'('s are re])eate(l until all colnnins and rows of B are ])ernintefl to (on'imi t a 
lowc'f triangnlar Idc^ek. T 1 k' ])ivot agenda algoritlini can l/e descril/ed mat liemat i< all\ 
as iollow": 

Procedure : /^\ o/-.\<!''//da 

Input : B. //. 

Output ; C. K. 
begin 

r ■— 1 ; // — ii: 

lot i I — 1. Ill dll J •— I ■— : L R — 0; 

lor I / — i.;/1 do 

lot I / = I.//1 do 

ili/i = III ilien ^7 —- / 1; I — I _{/); , » Init iali/.at ion • 

eiid lot 
end lot 

£ — -j - : 1 £ / < //. ^7 = (■.}; > Ifeeind the zero rows foi iceovt'r.x' » 

wliil' ' / £ //1 do 

CuLSiimh u. u. r.J. .C. K c 

Itiiw-Siunli II. II. r.l. J.C. R i; 

Buwji.Proci II. a. v. B. J. J. C. R. 3): 
end wliile 

end. 
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3. Column Singletons 


\\f ictci III C K ri' a jiaiiial |>i\iii am-inia il 1 here cxisian iih'h'x 1 :i /■ l " 
Mull tlia1 C, = Tv. = (I, lu'I I = •{ / ; /j 3; U. / Tv. I "}. 111! j ~ C. ll 

j J, ! = I. I hen rolunin / i" called a cohnnu ^nifih tan. IMacinu a culuiiin '-iiiijletun and 
llie 1 orre'.jKiiidini; luw Ui the liulii inn'-i una'-'iunei! ]M)''iiii>n ul ilie ])i\ui atienda li.e. 
C,: = Tv„ = 0 toy the larije^t ii) ie>ull' in aiipeiuiinu <ine addiliunal hiuei trifin[>ulai 
coluinn. Jn this iJiueeduie \\c rejx'ai the .search lor coluriiii sintileton-^ until all row 
and eulunin' are a'"'ii;n<‘d lo L'K’. oi nniil no coluinn siiuileioH'- are jound. 1 ho 
■'eai i ll I'll culuiiin ^iii'ileton^ can in- re[)ie-,enied mat heinai u a!!', a- lo!!o\' : 
Procedure : (' o/.s/na/ 

Input : //.II. r. J. ^7. 

Output : ii.I.J.C.K. 
begin 

whilei \ ZJ J — \ \ 7. " ' ■' 

I. - I. M ~17J: 

J — . n •— i; — ! : 

end while 

end. 

4. Row Singletons 

Let C and Tv Ire any partial pivot agenda, and let J, — {/ ; /i, (j.j ^ C. 1 < 
J < n}. For I ^ Tv. if |^7, j = 1. then row / is called a vow ■'■///(ilf Ion. Placing a row 
singleton and the corres)Jonding column to the left most unassigned position of the 
pivot agenda (i.e. C = Tv, = 0 for thesniallest r) results in aiipending one additional 
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Iilwci 11 iaii'iiil.ii Ill llii- |)r"<ic<lini- wc M-iii-al iln- Ini iii\', !■ a:' mi', 

all i<)\\ aiiii i iiliiiiiii- arc a-'''i‘;iii-(! i" c’' K >. <>1 iinii! im idvv 'IiilI( !iiIi' nM Imjia.' I in 
scarcli i'ur row siiiuicloiii' can l>c rcpn-'cnicd matlicinalicali\ as i'ollnw: 

Procedure ; 

Input : n. 11. r. I. J. 

Output : C, I. J . C. K. 
begin 

wliilcl I m : !i.7 = 1. \ S ' '' $ " '‘I" 

K - r. C - J ^ J- 
J. - J. 1,/! 7/ £ J; 

I ~ c - c - 1: 

(■ml wliilc 
end. 

('onsidcr t Im mat rix illusi rated in I i<>ur<'2. Ai»i)l\inu procedure r’o/_.S'in(;/to ild- 
matri.x re.sulis in placinti colntiin' ! ami 2 in the last two column' ol the ]iermuiei| 
matrix. .\ppl\iim procedure //ou.Si;,e'/ i,, the permuted matiix re'uli' iii placin'^ 
roW' (i. !l. II. ami 12 in the hist lour toW' of the permuted matrix. 1 he jx’rmnleo 
matrix i' illu'ttated in I'iuiiri 

5. The Bump 

.\ftpr applying tlie f 'oL.Singf and i^oux.Sin"/jrrocedures. the remaining rows. \i : 
1 < t < n. i ^ 7v}. and coluitms. {,/ ; J S 7 < ».j ^ C}. will eithei have jjZ,; > ] 
and jJjj > 1. or = U. For this section, we assume tliat |27di27,| ^ 0. These 

remaining rows and columns will form a nontriangular section of the permuted matrix, 
called the buiiijt. 
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COLU 

ROW IJl R 
1 1 
2 1 

3 le 

4 18 

5 16 

6 1 

7 9 


1 2 3 4 5 6 7 


9 10 11 12 13 14 15 16 17 16 19 2C 


e 

9 

10 

11 

12 

13 

14 

15 

16 
17 
16 
19 
2C 

COLP 


3 

1 

3 

1 

1 

3 

5 

e 

2 

9 

3 
7 

4 

I I I 
C 


X 

X 

X 


X 

X 

X 


X 

X 

X 


0 0 0 0 


9 10 11 12 13 14 15 16 17 16 19 20 
747 10 12 $666556 

00000 0 00000 c 


Figure 2: \ U* Kool^ aii manix. 

\W a pcniiui at ii^ii ol' tlu' Immp row- aii<l rohmiii'' that ii< yicM- uii)\ a 

jow '-mall i>lo(k.' that at»ovo i!h diaiional. iiit that niiiiimi/o- thf numiK': oi 

tioii/oi (M’lil ri< *- that appeal al >u\'e t h* • oliaiional )ij i he><- hKxk'. aiH 1 ' ii i > t hat < i juImk ■- 
t lie-Ne <’iii lie" 1 1 > a lew ewlinnU" that <*:\l <')m 1 a!x a e the (lia<juiial. < cilloi 1 i h* >////e 1 1 > 

ol>iaiii thi" juTimiiat iwii, -eleet a euhunii oi a >ei ol eohuim". whleh will iiit rothii «■ 
iIk' maxiiiiuin minihei ot row sinuletoii" wlien itMiiporarilv reni(>v(‘fi Iroiii the i>nnip. 
Let = !{/ : { £ J,. 0 < IvZ; ^ which i> called tije order fallt/ fun( tion wl 
column j. For a given k\ determine //.- for every column in the bumi>. If the maximum 
is greater than one. select the corresponding column for temporary removal using 
the maximum |JJ to break tie^. If the maximum /*,. equals one. only one row will lie 
affected i)y the temjiorary removal of this column. For this cas(‘. increase F to the 
minimum \J,\ greater than k and repeat the process. Let S denote the sequence of 
columns which are temporarily removed, (spikes). A column selection meiliod using 
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13 14 IL 16 17 18 IS 


11 0 11 
12 0 12 



X X X X X X X 
X X X X X X X 
X X X X X X X 
X XXX 


9 11 12 

c- 0 0 

9 11 12 


13 14 15 16 17 18 19 

12 5 6 6 6 5 5 

0 0 0 0 0 C 0 

2 0 0 1 0 1 0 


Figure 3 : Tli*' mairix al^T ii|»pi\iiiu tli<- ('(/LSiiii>l and Uiak ^Siaul proc^-diii^' 


itir Utily fiiiicntni can Kc (U'-^crilK'd luai iiciiiai icallx' a-- follow- 
Procedure : .\/a.\_7a//r 


Input : /. . n.I.C.S. 

Output : / . ' 
begin 

A ^ : 1 S ./ S .! 'z ./ r m \ \ 

wliile 1/// = ] vV A ^ O) cl<j 

for [j € A') do *- |{/ : / € I,. 0 < \J:\ < k}\\ 
IV *- max{/;,.(.y) : j € A'}: 
if( in = 1 ItlK'ii 

A' {.y : 1 < ,/ < n. tk(j) = 1}; 

<— n]in{|\7,i ; /• < \ J,\. 1 < ? < //}• 

end if 
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• •ml wliil'- 

— aimiiax-l I : / ^ A . /.' / ' = n, }; 

end. 


Block 



Figure 4 : k |)ariiiioti of a l)iniii). 


It flic I'uriip i- intf< a snial) ninnher of blocks rcjniioftod liy 

lowoi triaimtilar foiii|»oiioiit,\ii o.\aiiii>|o c>f t In'- sliurluro iliu^lraiod in Fiinif i. 
WitJjiij a bl(j(k. only tlio >jnk<','' aro j>orjiiin<'<! to haxo a zeio on tlio diasional. J lie 
bum]) procesjsing techniquo uses tlie MHx.TaUy procedure for .spike select ion. It liegins 
by setting k to the niinimum \J,\. A column is selected, using the Max.TaHy ])roce- 
dure. removed from the bump, and assigned to the spike sequence S. 1 his process is 
repeated until at least one row singleton is obtained. If the row singleton is unique, 
the corresponding pivot is placed in C,.(7v, ). Otherwise k is set to 1. a column. .-. 
is then selected using the Max.Tally |)rocedtire and placed in C.. Then one of the 
row singletons that lias a nonzero entry in column is selected, and placed in K . If 
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ij — / 11 i • 1. i.c. 1 III' nuinlii'i It! I'ow ''iii<;lci<iii-- iiav'iiit: iiuti/i'in <■ 111 : 11 '' iii <</i!i!i;i. - 
i' lliaii itiir. Wf lli< uiiil y Ut ])la('<' ' </ — 1 ' -'piki- iiiio lla |)!\i)' 

agrixla. Including all llic spiko'' into the pivot agenda isolates a hloi k. and de( i)ni- 
jxtses the Ituniji into smaller liloiks. 1 he piores'' is K'jx'aierl lor the idem ili< at ion <<! 
each Ithick. If B ha" a row of all zero', or a column of all zero', i.e. !J J = 0. 
then the de'ired lowei trianuula! 'tructiire i' not achi('\al)le. II thi' i' di'cox'ered. 
a reiuverv pioeedure. whi< h is d<''criliei| in the ne.vt section, will he inroked. 1 he 
hum]) jU'oressing can he descrihefl math<Mnati<ally as follows: 

Procedure ; Huin/j.Proc 
Input : I/. It. r. B. I. J. 3. 

Output : It. f.'B.J. 3 -C ■ ■ 

begin 

I — ( 1:5 — 

/■' — min { ■ J. — C’. \ < i 3 

if I I, ; ^7 = c/. 1 £/<//) = <:■' tV ii = <:> I t hen Reror'er.J \ k. n. 11 .'B.I .C .K. 3 r. 

while I (/,' > liit > (I I .V r < til (h< 

while I/.' > 1 I do ’ tem|)orarilv rcunove a set of column' - 

Mtix.'Iallyi /.. It. J. C. S. ■'!. 

I ~ I ~ 1: S. — J, — J, h'} v/.' -E Jo J. — o: 

I. — min{ ,J, : J. ^ o. \ < I < It \ : 

end while 

3^^ {1 : \J,\ =1.1 < I < //}: / » record the row singletons * j 

while [T' ^ o k ( > 0) do 

if] iJFj = 1 )then / ♦ a unique row singleton * / 

C I ^ J £ J'R, *— i € 3’’- 

Ji: - JMj} € J,: I, - o: c ^ c + 1: 

else 
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/• ^ 1 ; .\l II .1 -1 <111 i/i I.. I!. I .C. S. ' 

<1 \i ~ I ^ i (■ : ' ('(iH i 11M11 ■ / j!i • 

f K ^ , q ,F- J.; 

V.7, ■ - ^’7,, {.'} V /•' c ■ X- — O. I -— /' -f- 1; X" ■— .X . ■{ / ]; 

\vllll«' I 1 / > (I iV (/I tit) 

V — (/ - 1: 

iii (j > (I hIk'ii , » iiirliKic a s|)ikc into tho ag,onrla » 

C, 1 — S,: Tv, — ! F T. '■J,‘ — 1): 

T — J- •( /}: / — ' — 1: r — /• — I: 

('Uc ii ! t > 0 >V £ = o I t lini 

l{i‘i <>\ t'! Jji I/. II. B. J. C. Tv. £- 3 K 

I'lKl li 

cinl vTiilc ■ > (/ > (I * 

(‘inl il »■ [F. = 1 » 

3 — | / . J = I . 1 

I Mil I while ^ .F — Cl • 

/. — mill j ^'7 ; ^7 = o. 1 _ ' £ " 

it']' : 3 = o. I ' £ " ! = o y 3 = oi iIk’ii Ri'('t>\y'r3i h. ii .B.J .C .K. S. 3 
end while 


end. 


C'oM'iider the permuted malri.'; iDn-irated in Figure :F In tlii'- malrix. tlie una'-- 
signed rows* and columns- form the buni|). f sing the Max-Talh procedure, column 
]3 wa.« selected for temporary removal from the bump. This removal introduced two 
row singletons, rows s and 18. ( sing the Max.Tally procedure, column -f was selected 
to be assigned to Ce. and row 8 was selected to be assigned to TZ,,. This assignment 
introduced a new row singleton, row ];3. Similarly, rows 20 and 1'3 were assigned 
to Tv- and Tv^. where the <orresjjondine columns and the order of placement were 
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COL# 6 9 11 12 3 4 5 i 6 1C 14 15 16 1^ 18 19 20 13 I 2 

ROW IJ! R 

6 0 6 X 

9 0 9 X 

11 0 11 X 

12 0 12 X 

3 13 0 X X X X 

4 13 0 X X X X 

5 13 0 X X X X 

7 6 0 X X 

BIO X 

10 3 0 

13 2 0 

14 3 0 X 

15 5 0 XX 

16 2 0 

17 4 0 X X X X 

16 1 0 X 

19 3 0 X X X 

20 2 0 X 

10 1 

2 0 2 X 

COL# 6 9 11 12 3 4 5 7 8 10 14 15 16 17 18 19 20 13 1 2 

II! 0000514554566655800 0 

C 6 9 11 12 0 0 C C C C 0 0 0 C 0 0 0 C 1 2 

Kal TK OICOOOO 0 00100 

K« 2 TF 2 1 

Aft«r including cclimn 4 and row 8 into the agenda. 

K*1 TK 0 00000000101 

1 4 

rev «. t.r.* •9*r.i4 

C OCC0102010 


Figure 5: ilii' lnjiii|i ali<'! I'Miioviim f'f)huij)i !•< H' a '[lik''. 


(ificrDiincfl U'ijiu iIm- tally i'lUKtiuii. J lioc arc j!li)>i i at nl in I'iiiun 

Since = '1. li.c. the niiuilici ol row >iiiglcton> tha' lia\'c nu7i/ciu cntiic^ iii 

colninii l(i i^ '-’l. one s])ikc can i>c placc<l into the ])i\ol agciKla. Since ct^himn 1-5 
wa- tlie only spike, includins it into tin* ]>i\ol agenda deconiposerl the ljuinjt into twu 
small blocks, as illustrated in Figiin' (>. Since mininiuni \J^\ — 1 for the remaining 
rows, the Row.Singl procedure is invoked to assign these row singletons to the pivot 
agenda. The resulting matri.x is shown in Figure 7. Since minimum |j7i | is now greatf-r 
than one. the Bunip.Froc procedure was invoked again, and the resnlting matrix i-- 
illustrated in Figure S. By processing the bump, the number of nonzero entries above 
the diagonal was reduced to only fifteen entrie.s which are confined to font spike'- in 
two block'. 
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Figure 6: 1 li< iiiiiinx ali<'i (Iccoiiiixisinu a bKx k. 

6. Row Recovery 


I lie ali<i\c al'ii irit inii' ai<' <lc>iuii<''l i<> |)i'0(linc a sparse iaduii/at ion ol llic 
iiixci'i' i<| a linear |)iuuraiiiiiiiin; l>a'i'. Due to liiiite ])ie(i';)()n ul inaeliiiie'. luuiiii- 
ofl <'nui' oe( 111 will! h iiiav iea<l lu a "ci ul l^asii <(>luinn'- liax iiin a rov. ul zeiD'. 
When tlii> (X(ui>. ilu' ela''>i< llelleniian-Hariek algorilhin faiK, in llii^ "ceiiuii. 
we luesenl ilie recoverv algoiiihnn' whicii ensui<- tiial the inescriiied lower triangulai 
form is altainairle. while iiitroducine a replacement column having a single nonzeny 
entry. The ingeriuit\ is in the selection of the basic column to l>e replaced. 

Let 3 = {/ ; 1 < / < II. = o}. Let e, denote tlie column of an identitx 
matrix. Based on the structure of the original matrix, recovery can be done at three 
different stages of the algorithm. 

The first case for recovery occurs when there exists at least one zero row. and all 
the remaininc rows have lieen assigned. Since onlv n ~ 3] yrivots have been a'^'-igned. 
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16 
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18 
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0 

0 

0 

0 

0 

0 

1 



11 12 
0 0 
11 12 


4 20 16 13 14 16 7 3 5 

OCOOOOC54 
4 20 16 13 14 16 7 0 0 


8 10 15 17 19 

5 4 6 6 5 

C 0 0 0 0 


Figure 7: '] li<' niati ix afU'r ajtpiyinu ih#* Kou .SinsJ procodurr fov ili#' s^'ccmd liui**. 


!J ; = 0 for tli»- KMiiainiiia ajlunuis. and it is rlc^ar that tlinsc cohinnis shuuld \k' 
rri)la<(*d. frciiii tlii> cax* (an !><• matlioiuatirally dt^scriljrd a^ follow": 

Procedure: 

Input : n. (/. B. I. 2, 

# Output ; //. r. K, 2. 

begin 

C,. — 7 -E : 1 < /-■ < //. /■■ £ C,I, = o); 

^ Tvi, / 'E 2: 2 ^ 2 ,{/}: i/ — n — 1: A’ — 0: 

replace column j of B b\' e,: 

end. 


This case is illustrated in Figure 9. 

The second case for recovery occurs when there exists a zero row in the middle of 
processing the bump. The matrices permuted by the P 3 algorithm have the interesting 
property of s]>ikes ap])earing as properly nested sets, for more detaib see Arioli d 
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Figure 8: 1 he mairix alter )>n)((‘>si!m tlie la>i hum]). 

(il. r. .\ 'pike liH' liccii <lisoi\<>ic(l for tlii" r<“<<)\'i'i\' case, wliidi (auropumi- lo tiu' 
outer iiio'l nc'lfd x t. aixt t iiat wiil never lx- placed into i he jtivot a»enda. 1 hi" >piki' 
i" al\va\> re])laeed hy the appropriate coltuiiii of th<‘ idenlil\' uiaiiix. l{(xuver> Iidiu 
the ‘^eeoiifl ease can l)e ntatheinatiealh de>eri!)ed as iollows: 

Procedure ; Hecovfi -2 
Input : II. ii.'B.J.S. Z. 

Output : it. B.c’. 7\. 5. 3. 
begin 

C,i *— J Sj: TZu ^ i € Z-: 

S^S\S,: Z ^ Z\{t}: f/-</-]: 

replace column j of B by e,; 

end 

This case is illustrated in FiEttre 10. 
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Figure 9: An ''Xiiinplc <>! rccovcrN ( ii'c 1. 

l iii' Uiird I'n'i' l<<i it(u\'''r\ .xcui^ wlim thcic cxi-l' n /cru row ai tlu' cinl ut 
])M>< r"'iiix I lie l)Uiii|). i'ni tl,;' < a-*', t in't'' a'r< . 'Z- — S uiiii^'iuiii'il ( i)liiiiiii' ainl a 
hfi" Itccii (li>f'>\ciril iliai will never he j)lared inU) ilie acenda. Similar in 

tlie ''(’(■(/ikI ea<e, t l)i> s|>ike is always reulared l.y i Ik- ajjjiropriaie rolmim of i lie ifleii’ iiy 
# iiialri.x iisiiiu [Moeediire U< <u\<'rj2. Kvemiially. w<' will roach a siiiiaiion similar 

the lirsl case, where jirucedure Rrunri^J will he invoked, Tlii" case i- ilhislia1 ed in 
Figure 1 1. 


7. Summary and Conclusions 


We believe that the ideas presented in this manuscript close a gap in the literature 
concerning the preassigned pivot agenda algorithms. .All of these algorithms assume 
an input matri.x which is nonsingular, that may not he the case every time thes<' 
procedure^ are called. This manus<-ri];t presents a \ariant cd which icsult-- in a 
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Figure 10; An pxaiii]))e of rerovory ( ^>0 2. 


relialile and eftifient iiietliod to recover from all types of numerical instabilities. 1 he^e 
ifleas ran be easily adajited for tin* aleoritlmis of [1], [3j. [Ij. [f>j. and [7j. 
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